# devtools::install_git('https://git.oschina.net/seifer_08ms/ROzone2.git') 
# 中文windows下R安装在c盘program file下此句会报错，
#建议R改成安装在d盘非中文无空格路径
library(ROzone2)
install.gdal() ## 默认参数兼容性更好，默认安装在c盘gdalwin下面。删除此文件夹会重新下载gdal二进制工具。
library(raster)
library(dplyr)
library(stringr)
# 数据放在当前路径的data文件夹下面，后缀名是nc
data.dsn.list<-list.files(path = 'data',pattern = '.nc$',recursive = T,full.names = T,include.dirs = T)
print(data.dsn.list)# 确认所有nc文件
#根据输入数据生成输出数据路径，把nc替换为tif
dst.dsn.list<-str_replace(data.dsn.list,'\\.nc$','\\.tif')
print(dst.dsn.list)# 确认nc文件保存后的路径
# ras.stacks<-stack(data.dsn.list%>%head)
# 实际调用的命令格式
# gdalwarp -q -cutline /tmp/wcl/chn.shp -crop_to_cutline -of GTiff /tmp/wcl/nc1998-2014-allcomposition/GlobalGWR_PM25_GL_199801_199812-RH35.nc /tmp/res.tif

library(gdalUtils)
mask.dsn<-'data/chn.shp'# 中国边界shp
library(rgdal)
# chn.shp<-readOGR(mask.dsn,layer = ogrListLayers(mask.dsn)[1],stringsAsFactors = F)
mask.prj4<-proj4string(china.pop)

input.nc.prj4<-raster(data.dsn.list[1])%>%crs
china.pop.proj<-china.pop%>%spTransform(input.nc.prj4)
# ogrInfo(mask.dsn)$p4s # 
##注意 不是ogrinfo命令.ogrinfo是gdalUtility的命令，ogrInfo是rgdal命令
##############################################################################################
# zonal statistics
zonal.shp<-vector(mode = 'list',length(data.dsn.list))
for (i in 1:length(data.dsn.list)){ # 循环调用gdal工具的命令切割
    gdalwarp(srcfile = data.dsn.list[i],cutline = mask.dsn,of='GTiff',tr=c(0.08,0.08),
             s_srs =input.nc.prj4, t_srs=input.nc.prj4,dstnodata = -9999,
             crop_to_cutline = T,dstfile = dst.dsn.list[i],verbose = T,stats = T,overwrite = T)
    
}
tif.ras<-raster(dst.dsn.list[1])
res.test<- res(tif.ras)[1]
zonal.ras<-rasterize.gdal(china.pop.proj[,'county_ID'],
                          'county_ID',
                          res=res.test,extent = extent(tif.ras))
#############################################################################################
# merge all data into single table 
o3.stack.3<-lapply(dst.dsn.list,raster)
names(o3.stack.3)<-str_replace(data.dsn.list,'\\.nc$','')%>%basename
zonal.res.all<-lapply(X=o3.stack.3,FUN = zonalStat,
                      zonal.ras=zonal.ras,prefix='z')%>%
    do.call('cbind',.)%>%
    mutate(z=.[[1]])%>%
    dplyr::select(-ends_with(".z"))%>%
    dplyr::select(z,everything())
###############################################################################################
## compute mean value of all layers
zonal.res.all$allmean<-
    apply(zonal.res.all%>%
              dplyr::select(ends_with('mean')),
          MARGIN = 1,
          FUN = mean,
          na.rm=T)
#  join table
china.pop.all<-china.pop.proj
china.pop.all@data<-china.pop@data%>%
    left_join(zonal.res.all,
              by = c("county_ID"="z")  )

china.pop.plot<-china.pop.all
# drop some columns which names are too long to be handled with ESRI Shapefile
selectnames<-names(china.pop.all)[
    (china.pop.all%>%names%>%length-length(data.dsn.list)):(names(china.pop.all)%>%length)]
china.pop.plot@data<-china.pop.all@data%>%dplyr::select(one_of(c('county_ID',selectnames)))
names(china.pop.plot)[2:(names(china.pop.plot)%>%length-1)]<-paste0('val',1:length(dst.dsn.list))
writeOGR(china.pop.plot,'data/chinapop.shp',driver = 'ESRI Shapefile',layer='chinapop',
         overwrite_layer = T,delete_dsn = T)
##################################################################################################
library(tmap)
qtm(china.pop.plot,fill='allmean',title = 'Distribution of PM2.5',
    fill.title='Legend\n(μg/m3)', style="gray")
